home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Pascal Super Library
/
Pascal Super Library (CW International)(1997).bin
/
MATH
/
PIBSIGS
/
GAMMAIN.PAS
< prev
next >
Wrap
Pascal/Delphi Source File
|
1986-06-22
|
6KB
|
189 lines
(*-------------------------------------------------------------------------*)
(* GammaIn -- Incomplete Gamma integral *)
(*-------------------------------------------------------------------------*)
FUNCTION GammaIn( Y, P : REAL;
Dprec : INTEGER;
MaxIter : INTEGER;
VAR Cprec : REAL;
VAR Iter : INTEGER;
VAR Ifault : INTEGER ) : REAL;
(*-------------------------------------------------------------------------*)
(* *)
(* Function: GammaIn *)
(* *)
(* Purpose: Evaluates Incomplete Gamma integral *)
(* *)
(* Calling Sequence: *)
(* *)
(* P := GammaIn( Y, P : REAL; *)
(* Dprec : INTEGER; *)
(* MaxIter : INTEGER; *)
(* VAR Cprec : REAL; *)
(* VAR Iter : INTEGER; *)
(* VAR Ifault : INTEGER ) : REAL; *)
(* *)
(* Y --- Gamma distrib. value *)
(* P --- Degrees of freedom *)
(* Ifault --- error indicator *)
(* *)
(* P --- Resultant probability *)
(* *)
(* Calls: *)
(* *)
(* None *)
(* *)
(* Remarks: *)
(* *)
(* Either an infinite series summation or a continued fraction *)
(* approximation is used, depending upon the argument range. *)
(* *)
(*-------------------------------------------------------------------------*)
CONST
Oflo = 1.0E+37;
MinExp = -87.0;
VAR
F: REAL;
C: REAL;
A: REAL;
B: REAL;
Term: REAL;
Pn: ARRAY[1..6] OF REAL;
Gin: REAL;
An: REAL;
Rn: REAL;
Dif: REAL;
Eps: REAL;
Done: BOOLEAN;
LABEL 9000;
BEGIN (* GammaIn *)
(* Check arguments *)
Ifault := 1;
GammaIn := 1.0;
IF( ( Y <= 0.0 ) OR ( P <= 0.0 ) ) THEN GOTO 9000;
(* Check value of F *)
Ifault := 0;
F := P * LN( Y ) - ALGama( P + 1.0 ) - Y;
IF ( F < MinExp ) THEN GOTO 9000;
F := EXP( F );
IF( F = 0.0 ) THEN GOTO 9000;
(* Set precision *)
IF Dprec > MaxPrec THEN
Dprec := MaxPrec
ELSE IF Dprec <= 0 THEN
Dprec := 1;
Cprec := Dprec;
Eps := PowTen( -Dprec );
(* Choose infinite series or *)
(* continued fraction. *)
IF( ( Y > 1.0 ) AND ( Y >= P ) ) THEN
BEGIN (* Continued Fraction *)
A := 1.0 - P;
B := A + Y + 1.0;
Term := 0.0;
Pn[ 1 ] := 1.0;
Pn[ 2 ] := Y;
Pn[ 3 ] := Y + 1.0;
Pn[ 4 ] := Y * B;
Gin := Pn[ 3 ] / Pn[ 4 ];
Done := FALSE;
Iter := 0;
REPEAT
Iter := Iter + 1;
A := A + 1.0;
B := B + 2.0;
Term := Term + 1.0;
An := A * Term;
Pn[ 5 ] := B * Pn[ 3 ] - An * Pn[ 1 ];
Pn[ 6 ] := B * Pn[ 4 ] - An * Pn[ 2 ];
IF( Pn[ 6 ] <> 0.0 ) THEN
BEGIN
Rn := Pn[ 5 ] / Pn[ 6 ];
Dif := ABS( Gin - Rn );
IF( Dif <= Eps ) THEN
IF( Dif <= ( Eps * Rn ) ) THEN Done := TRUE;
Gin := Rn;
END;
Pn[ 1 ] := Pn[ 3 ];
Pn[ 2 ] := Pn[ 4 ];
Pn[ 3 ] := Pn[ 5 ];
Pn[ 4 ] := Pn[ 6 ];
IF( ABS( Pn[ 5 ] ) >= Oflo ) THEN
BEGIN
Pn[ 1 ] := Pn[ 1 ] / Oflo;
Pn[ 2 ] := Pn[ 2 ] / Oflo;
Pn[ 3 ] := Pn[ 3 ] / Oflo;
Pn[ 4 ] := Pn[ 4 ] / Oflo;
END;
UNTIL ( Iter > MaxIter ) OR Done;
Gin := 1.0 - ( F * Gin * P );
GammaIn := Gin;
(* Calculate precision of result *)
IF Dif <> 0.0 THEN
Cprec := -LogTen( Dif )
ELSE
Cprec := MaxPrec;
END (* Continued Fraction *)
ELSE
BEGIN (* Infinite series *)
Iter := 0;
Term := 1.0;
C := 1.0;
A := P;
Done := FALSE;
REPEAT
A := A + 1.0;
Term := Term * Y / A;
C := C + Term;
Iter := Iter + 1;
UNTIL ( Iter > MaxIter ) OR ( ( Term / C ) <= Eps );
GammaIn := C * F;
(* Calculate precision of result *)
Cprec := -LogTen( Term / C );
END (* Infinite series *);
9000: (* Error exit *)
END (* GammaIn *);